/**
 * @file   Interp.h
 * @author xingyifan
 * @date   2021/10/25 10:56:17
 *
 * @brief  struct InterpCondition, class Polynomial, class NewtonInterp
 *
 *
 */

#include "Math.h"

using namespace std;

#ifndef __INTERP__
#define __INTERP__

typedef struct InterpCondition
{
    vector<double> points;                   //插值点
    map<double, vector<double>> interpSites; //插值点对应函数值
} INTERPC;

class Polynomial
{
private:
    vector<double> coefficient; //系数向量

public:
    Polynomial() { coefficient.push_back(0); }
    Polynomial(vector<double> _coefficient) : coefficient(_coefficient) {}

    //获得函数值
    double caculate(double);

    Polynomial derivation();
    Polynomial derivation(u_int);       //求导
    Polynomial integration();           //不定积分
    double integration(double, double); //定积分

    static Polynomial Chebyshev(u_int); //获取Chebyshev多项式

    Polynomial operator+(const Polynomial &poly)
    {
        Polynomial res = *this;
        int index = 0;
        int this_size = get_degree() + 1;
        int that_size = poly.get_degree() + 1;

        while (index < MIN(this_size, that_size))
        {
            res.coefficient[index] += poly.coefficient[index];
            index++;
        }

        while (index < that_size)
        {
            res.coefficient.push_back(poly.coefficient[index]);
            index++;
        }
        return res;
    }

    void operator+=(const Polynomial &poly)
    {
        *this = *this + poly;
    }

    Polynomial operator-(const Polynomial &poly)
    {
        Polynomial res = *this;
        int index = 0;
        int this_size = get_degree() + 1;
        int that_size = poly.get_degree() + 1;

        while (index < MIN(this_size, that_size))
        {
            res.coefficient[index] -= poly.coefficient[index];
            index++;
        }

        while (index < that_size)
        {
            res.coefficient.push_back(-1 * poly.coefficient[index]);
            index++;
        }
        return res;
    }

    void operator-=(const Polynomial &poly)
    {
        *this = *this - poly;
    }

    Polynomial operator*(double k)
    {
        Polynomial poly = *this;
        for (int i = 0; i < (int)get_degree() + 1; i++)
        {
            poly.coefficient[i] *= k;
        }
        return poly;
    }

    Polynomial operator*(const Polynomial &poly)
    {
        Polynomial res;
        int m = get_degree() + 1;
        int n = poly.get_degree() + 1;
        res.coefficient.resize(m + n);
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n; j++)
            {
                res.coefficient[i + j] += coefficient[i] * poly.coefficient[j];
            }
        }
        // clean zero
        reverse(res.coefficient.begin(), res.coefficient.end());
        vector<double>::iterator it = res.coefficient.begin();
        while (it != res.coefficient.end())
        {
            if (*it == 0)
            {
                it = res.coefficient.erase(it);
            }
            else
            {
                break;
            }
        }
        reverse(res.coefficient.begin(), res.coefficient.end());
        return res;
    }

    void operator*=(double k)
    {
        *this = *this * k;
    }

    void operator*=(const Polynomial &poly)
    {
        *this = *this * poly;
    }

    void operator=(const Polynomial &poly)
    {
        coefficient = poly.coefficient;
    }

    bool operator==(const Polynomial &poly)
    {
        if (coefficient == poly.coefficient)
        {
            return true;
        }
        return false;
    }

    Polynomial operator/(double k)
    {
        Polynomial poly = *this;
        for (int i = 0; i < (int)get_degree() + 1; i++)
        {
            poly.coefficient[i] /= k;
        }
        return poly;
    }

    friend ostream &operator<<(ostream &out, const Polynomial &poly)
    {
        out << "p(x) = " << poly.coefficient[0];

        for (int i = 1; i < (int)poly.coefficient.size(); i++)
        {
            out << " + " << poly.coefficient[i] << "x^" << i;
        }
        out << endl;
        return out;
    }

    //获取多项式的阶
    u_int get_degree() const
    {
        return coefficient.size() - 1;
    }

    //获取系数
    vector<double> get_coefficient() const
    {
        return coefficient;
    }
};

class NewtonInterp
{
private:
    INTERPC interpC_;                            //插值条件
    Polynomial interPoly_;                       //插值多项式
    vector<vector<double>> tableOfDividedDiffs_; //差分表

    void caculateTable(int, int);                                            //计算差分表
    void caculatePoly();                                                     //由差分表计算插值多项式
    static Polynomial recursionPoly(double, double, Polynomial, Polynomial); //递推多项式
    static Polynomial Neville_Aitken(INTERPC &, int, int);                   //内部递归的算法过程

    static double recursion(double, double, double, double, double); //递归计算过程
    static double Neville_Aitken(double, INTERPC &, int, int);       //递归计算x处的值

public:
    NewtonInterp() {}

    void overwriteInterp(INTERPC); //重新插值
    void incrementInterp(INTERPC); //递推插值

    static double Neville_Aitken(double, INTERPC &); // Neville Aitken算法
    static Polynomial Neville_Aitken(INTERPC &);     // Neville Aitken算法
    void Chebyshev(u_int, double (*)(double));       // Chebyshev插值

    friend ostream &operator<<(ostream &out, const NewtonInterp &interp)
    {
        out << interp.interPoly_ << endl;
        for (int i = 0; i < (int)interp.get_size(); i++)
        {
            vector<double> arr(interp.tableOfDividedDiffs_[i]);
            out << "{";
            for (int j = 0; j < (int)arr.size(); j++)
            {
                out << arr[j] << ", ";
            }
            out << "\b\b}" << endl;
        }
        return out;
    }

    u_int get_size() const //插值点的数量
    {
        return tableOfDividedDiffs_.size();
    }

    INTERPC get_condition() const //插值条件
    {
        return interpC_;
    }

    Polynomial get_poly() const //获取插值多项式
    {
        return interPoly_;
    }
};

#endif
